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ABSTRACT 



Aims. We investigate the effect of including a proper energy balance on the interaction of a low-mass planet with a protoplanetary 
disk. 

Methods. We use a three-dimensional version of the RODEO method to perform hydrodynamical simulations including the energy 
equation. Radiation is included in the flux-limited diffusion approach. 

Results. The sign of the torque is sensitive to the ability of the disk to radiate away the energy generated in the immediate surroundings 
of the planet. In the case of high opacity, corresponding to the dense inner regions of protoplanetary disks, migration is directed 
outward, instead of the usual inward migration that was found in locally isothermal disks. For low values of the opacity we recover 
inward migration and show that torques originating in the coorbital region are responsible for the change in migration direction. 
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1. Introduction 

Migration has become a standard ingredient in planet formation 
theory, providing an explanation for the extrasolar giant planets 
that were found orbiting very close to their central star, the so- 
called Hot Jupiters. It is believed that these planets were formed 
at a distance of several Astronomical Units (AU), after which 
they moved in due to tidal interaction with the surrounding pro- 
toplanetary disk. 

Depending on the mass of the planet, three different migra- 
tion modes can be distinguished: 



1. 



3. 



Type I migration (G oldreich & Tremainell 9711 : IWardl 19971) . 
valid for low-mass planets. The disk response is linear in 
this case, and the resulting torques ca n be calculated in a 
semi-analytic way (Tanaka et al. 2002). Planets that are less 
massive than a few Earth masses are subject to this mode of 
migration. 

Type II migrati on, for which the pl anet orbits in a deep an- 
nular gap (Lin & Papaloizou 1986). In this case, the planet 
is tidally locked inside the gap, and the time scale for in- 
ward migration equals the viscous time scale, which approx- 
imately equals 10 6 yr. This mode of migration holds for plan- 
ets comparabl e to o r more massive than Jupiter. See, how- 
ever, |Rafiko3 {2002 ) on gap formation for low-mass planets. 
In between these two regimes, interesting things happen. 
The onset of non -linear behavior can influence the torque 
to a large extent (Masse t et alJl200 6a). However, when the 
disk is massive enough, the planet enters a new migration 
mode that can be very fast and that relies on dynamical 
corotation torques (Masset & Parjaloizou 2003; Artvmowicz 
l200l . This is called Type III migration, and it holds for plan- 
ets that open up a partial gap. 
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In this Letter, we focus on deeply embedded low-mass planets 
and, therefore, on the regime of Type I migration. 

The time scale for Type I migration is inversely p roportional 
to the mass of the disk and the planet dTanaka et alJ l2002). and 
it can be much shorter than the disk lifetime of approximately 
10 7 yr. A planet of 1 Earth mass (1 M e ), for example, em- 
bedded in the minimum mass solar nebula (MMSN) at 5 AU, 
would migrate into the central star within 10 5 yr. Within the 
widel y accepted cor e-accretion model of giant planet formation 
(PollacketaD[l996), the gas giants also pass through a stage 
where the core is a few M ffi , which thus should disappear into the 
central star even faster. Therefore, some sort of stopping mecha- 
nism is required for planets to survive Type I migration. 

Proposed stopping mech anisms include magnetic turbul ence 
(Nelson & Papaloizou 2004), a toroidal magnetic field (|Ter auem 
2003 ), and an inner hole in the disk ( Kuchner & Lecat 2002; 
Mas set et alJl2006a) . However, most analytical and numerical 
work on planet migration has focused on disks with a fixed ra- 
dial temperature profile. For these disks, hydrodynamical sim- 
ulations can repro duce Type I migration in both two and three 
dimensions dD' Angelo et all2002l2003l:lBate et all2003l) . 

The isothermal assumption is only valid if the disk can radi- 
ate away all excess energy efficiently. When the disk is optically 
thick, the radiative energy flux through the surface of a sphere 
of radius H around the planet is given by Fr = ctT 4 /t, where cr 
is the Stefan-Boltzmann constant and t is the optical depth over 
a distance H. Therefore, when the opacity increases, less energy 
can be radiated away. This makes the cooling time scale much 
longer than the dynamical time scale in the dense inner regions 
of protoplanetary disks, and therefore the isothermal assumption 
is likely to be invalid there. 

One may wonder what effects releasing the isothermal as- 
sumption may have on the torque. The effects of a more realistic 
temperature profile on the Lindblad torques were investigated 
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bv Uang-Condell & Sasselovl (120051) and iMenou & Good man 
(2004), who found that the migration rate is very sensitive to the 
temperature and opacity structure and that in some cases migra- 
tion may be very slow compared to the isothermal Type I case. In 
this Letter, we present the first results of including a detailed en- 
ergy balance on the migration behavior of low-mass planets in a 
radiation-hydrodynamical context. In Sect.[5]we discuss the disk 
model, in Sect.[3]we present the results, and we give a discussion 
of the results in Sect.0] We then conclude in Sect. [5] 

2. Disk model 

We work in spherical coordinates (r, 8, (p) with the central star 
in the origin, and the unit of distance is the planet's orbital 
radius. In the plots we also use a Cartesian coordinate frame 
with the central star in the origin and the planet located at 
(x,y,z) = (-1,0,0). The planet stays on a fixed circular orbit 
throughout the simulation. The disk is three-dimensional, and 
the computational domain is bounded by r = 0.4, r = 2.5, 
9 - n/2, and 9 = n/2 - 5h/2, where h is the relative scale 
height of the disk (h = H/r, with H the pressure scale height). 
The disk spans the full 2n in azimuth. For the radial and merid- 
ional^ boundari es, we use non-r eflective boundary conditions 
( Paardekooper & Mellema 2006), except for the boundary at 
9 = n/2, which is fully reflective. The azimuthal boundary con- 
ditions are periodic. This domain is covered by a grid with 256 
radial cells, 768 azimuthal cells, and 16 meridional cells. This 
configuration leads to cubic cells near the planet. On top of this 
main grid, we put 5 levels of adaptive mesh refinement (AMR) 
near the planet. Each level increases the local resolution by a 
factor of 2, which leads to a local resolution of Ar w 0.00025. 
We focus on a planet of 5 M ffl located at 5 AU from the cen- 
tral star, for which the Hill radius is then covered by 70 cells. 
While this resolution is not needed for accurate torque calcu- 
lations, it is neede for correctly modeling the planet's envelope 
and for the process of accretion, which we will consider in a 
future publication. In the torque calculations, we neglect the 
material orbiting inside the Hill sphere. For deeply embedded 
objects, material co uld be bound in the smaller Bondi sphere 
(Masset et al. 2006b), so that we may exclude a region that is 
too large. However, we checked that the results are not sensitive 
to this. 

The disk is taken to be inviscid, so as to filter out effects of 
viscous heating. The density follows a power law initially with 
index -3/2. We take h = 0.05, a typical value used in simu- 
lations of disk-planet interaction. This sets the temperature at 
approximately 50 K at 5 AU. The disk is in Keplerian rotation 
initially with a slight correction for the radial pressure gradient. 
The potential contains terms due to the star, the planet, and the 
acceleration of the coordinate frame due to the planet and the 
disk. The last two contributions are of negligible importance for 
a planet of 5 M e . We smooth the potential of the planet over 2 
grid cells. 

We solve the flow equatio ns using a three-d imensio nal ver - 
sion of the RODEO method (Paarde kooper & Mellemal l2006). 
with an energy equation included. Radiative transfe r is treated 
in the flux-limited diffusion approach ( Levermore & Pomraning 
119811) . together with the flux limiter of lKlevI dl989l) . The three- 
dimensional hydrodynamics solver (including AMR) was tested 
through isothermal disk-planet interaction, which showed that 
we could reproduce Type I and Type II migration in the rel- 
evant mass regime. The radiative transfer module was tested 
against multi-d imensional d iffusion problems. We used opacity 
data from lBell & Linl JI994). In the temperature range of inter- 
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Fig. 1. Total torque on a 5 M ffi planet as a function of time for 
three different midplane densities, together with the isothermal 
res ult. The torques are normalized to the analytical value found 
bv lTanaka et alJ (120021) . which is reproduced by the isothermal 
simulation. For high densities (and thereby for high opacities) 
the torque becomes positive, indicating outward migration. 

est, the opacity is proportional to the density and temperature 
squared. We vary the initial midplane density to study different 
cooling regimes while keeping the initial temperature fixed. 

The code is parallelized using MPI, and the simulations 
were run on 126 1.3 GHz processors at the Dutch National 
Supercomputer. Each radiation-hydrodynamical model took ap- 
proximately 4 days of computing time. 

3. Results 

In Fig.^ we show the evolution of the total torque on our 5 M ffi 
planet for three different midplane densities, where po = 10" 11 
g cm" 3 corresponds to the MMSN at 5 AU. For this high density, 
the torque is positive, indicating outward migration. The abso- 
lute value of the torque is approximately a factor 2 lower than 
the Type I torque. For lower densities, and therefore for lower 
opacities, we recover inward migration; and for the lowest value 
of the density, we approach inward Type I migration of isother- 
mal disks. 

To find out where this positive torque originates, we show 
the radial torque profile close to the planet in Fig. [2] An embed- 
ded planet excites two spiral waves into the disk. The inner spiral 
wave is responsible for the positive bump near r = 0.95, while 
the outer spiral wave creates the negative bump near r = 1.05. 
The location and the strength of these bumps varies between the 
different runs, but the total Lindblad torque remains negative for 
all simulations. For po = 10" 11 gem" 3 , the magnitude of the 
Lindblad torque is a factor of 2 smaller compared to the isother- 
mal case. This implies that it is the corotation region, located be- 
tween r = 0.95 and r = 1.05, that is responsible for the change 
in sign of the total torque. 

We can pin down the origin of the positive corotation torque 
further by looking at the vertical torque distribution. In Fig. [3] 
we compare an isothermal run to an adiabatic simulation. The 
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Fig. 2. Radial torque distribution for a 5 M e planet for two dif- 
ferent densities, together with the result for a locally isothermal 
equation of state. Both the inner wave and the outer wave are 
reduced in strength for higher densities, but the major difference 
comes from the corotation region. 
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Fig. 4. Temperature at z = 0.02, slightly above the Hill sphere 
of the planet. Overplotted is the direction of the velocity field, 
clearly showing the horseshoe orbits. The point where the two 
horseshoe legs meet is a region of compression, meaning the 
temperature is slightly higher at that location. 
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Fig. 3. Vertical torque distribution for an adiabatic simulation 
with no radiative cooling compared to the locally isothermal 
result. The difference starts high above the Hill sphere of the 
planet, which is approximately at z = 0.017. 



latter does not include radiative cooling, so the effect of ineffi- 
cient cooling is even more pronounced. The important thing to 
note is that the difference in Fig. [5] starts well above the Hill 
sphere of the planet, the edge of which is located at approxi- 
mately z = 0.017; therefore, we can filter out the effects of heat- 
ing deep within the planetary envelope and look at temperature 
differences above z = 0.017. 



We show a cut through z = 0.02 in Fig.@] showing the tem- 
perature and the direction of the velocities. The horseshoe orbits 
clearly stand out, and where the two horseshoe legs meet, the 
gas is compressed and the temperature rises. Due to the slightly 
sub-Keplerian velocity of the gas, the two horseshoe legs meet 
behind the planet (see Fig.|3}. As the disk tries to maintain pres- 
sure equilibrium, the density will be lower in this region. We 
verified that the temperature rise is enough to account for the 
positive torque by verifying that 

r GM p dp , 

i P -^r p x(r-r p )d 3 r, (1) 

with dp = -podT/To, is within the same order of magnitude as 
the total torque on the planet. Equation (1) reproduced the total 
torque to within 25 %, which shows that indeed the tempera- 
ture change is strong enough to induce a density asymmetry that 
is capable of producing the observed positive torque. It is this 
density asymmetry that leads to a large, positive contribution to 
the torque. Only if the disk can radiate away this heat efficiently 
does the total torque become negative again. 

4. Discussion 

The corotation regi on is very sensitive to local conditions 
( Mas sel|l200l[ 120021) . including the shape of the planetary en- 
velope. This may affect the contribution to the total torque. 
However, we have found that the positive corotation torque is 
a robust effect, showing no strong dependence on numerical res- 
olution, for example. Because we take a smoothing length of 2 
grid cells, we have also found no dependence on the exact form 
of the planetary potential. 

We estimate that the transition between inward and outward 
migration happens around 15 AU in the MMSN, based on the 



4 



Sijme-Jan Paardekooper and Garrelt Mellema: Halting Type I planet migration 




V 



Fig. 5. Temperature at x — -1 (or r — 1), showing the hot plume 
behind the planet (y > 0, the planet moves towards the left of the 
figure). This plume extends all the way to z = 0, but there the 
heating inside the planetary envelope dominates the temperature 
structure. For z > we show a simulation without a global pres- 
sure gradient (i. e. the gas orbits at the Kepler velocity), while 
for z < we show a simulation with a global radial pressure 
gradient (p oc pT oc r~ 2 ). The difference in the position of the 
warm plume is apparent for both simulations. The vertical line 
at y = is there to guide the eye. The color scale has been chosen 
to focus on the warm plumes. 



simulations with different opacities. Although we parametrized 
the opacity using the local gas density, it is the opacity that mat- 
ters. If the opacity is lowered, as when it is due to grain growth 
(Dullemond & Dominik 2005), the transition radius will be lo- 
cated farther inward. 

The computational costs of the radiation-hydrodynamical 
simulations is too high to run simulations for hundreds of orbits. 
Therefore, effects that happen on the librat ion time scale, such 
as saturation of the corotation torque (Ogilvie & Lubowl l2003). 
are not captured by these simulations. Adiabatic runs showed 
that indeed saturation occurs on the libration time scale, which 
considerably reduces the magnitude of the positive torque. This 
is another confirmation that is indeed the corotation region that 
is responsible for the positive torque. This does not mean, how- 
ever, that the positive torque should disappear on a libration time 
scale. When the radiation diffusion time scale is shorter than the 
libration time scale, the asymmetry in temperature, and therefore 
in density, can be maintained. In the MMSN, this is the case for 
approximately r > 1 AU. 

For lower-mass planets, the horseshoe region shrinks, so the 
positive torque should decrease in strength. However, in this case 
the Lindblad torques also decrease in magnitude. We have veri- 
fied that even for a planet of 0.5 M ffi the total torque is still posi- 
tive for po = 10 -11 g cm -3 . Further study is required to pin down 
the exact dependence of the torque on planetary mass. We expect 
that in the high-mass regime, when the planet starts to open up a 
gap and when the density around the planet decreases by several 



orders of magnitude, cooling will always be efficient. Therefore, 
Type II migration should not be very d ifferent from the isother- 
mal case (see also Klahr & Kiev 2006). 

5. Conclusions 

In this Letter, we have presented the first radiation- 
hydrodynamical simulations of low-mass planets embedded in 
a protoplanetary disk. We find that the migration behavior of a 
5 M ffl planet depends critically on the ability of the disk to radi- 
ate away the heat generated by compression close to the planet. 
When radiative cooling is not efficient, the torque on the planet 
is positive, indicating outward migration. Analysis of the torque 
distribution showed that it is the coorbital region that generates 
the large positive contribution to the torque, while the Lindblad 
torque is reduced in strength, but is still negative. Local, adia- 
batic simulations showed that compression at the turn-over point 
of the horseshoe orbits creates a warm region behind the planet. 
As the disk tries to maintain pressure equilibrium, the density 
behind the planet will become lower, which leads to a positive 
corotation torque. This happens in the dense inner parts of pro- 
toplanetary disks, where the opacity is high. In regions of low 
opacity, we recover inward Type I migration. Therefore, low- 
mass planets do not migrate inward all the way to the central star, 
bun instead stop migrating when radiative cooling becomes inef- 
ficient. In the MMSN, the transition occurs at approximately 15 
AU. As the disk slowly accretes onto the central star, this radius 
will move inward; as a result, low-mass planets will also even- 
tually continue their inward migration, but on the much longer 
viscous time scale. This way, low-mass planets are saved from 
fast inward Type I migration into the central star. 
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